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Abstract. We derive the mass model of the Milky Way (MW) using a cored dark matter (DM) halo profile and recent data. The method 
used consists in fitting a spherically symmetric model of the Galaxy with a Burkert DM halo profile to available data: MW terminal 
velocities in the region inside the solar circle, circular velocity as recently estimated from maser star forming regions at intermediate 
radii, and velocity dispersions of stellar halo tracers for the outermost Galactic region. The latter are reproduced by integrating the Jeans 
equation for every modeled mass distribution, and by allowing for diff'erent velocity anisotropics for diff'erent tracer populations. For 
comparison we also consider a Navarro-Frenk- White profile. We find that the cored profile is the preferred one, with a shallow central 
density of p// ~ 4 x lO^Mo/kpc^ and a large core radius Rh ~ lOkpc, as observed in external spirals and in agreement with the mass 
model underlying the Universal Rotation Curve of spirals. We describe also the derived model uncertainties, which are crucially driven 
by the poorly constrained velocity dispersion anisotropics of halo tracers. The emerging cored DM distribution has implications for 
the DM annihilation angular profile, which is much less boosted in the Galactic center direction with respect to the case of the standard 
ACDM, NFW profile. Using the derived uncertainties we discuss finally the limitations and prospects to discriminate between cored 
and cusped DM profile with a possible observed diff'use DM annihilation signal. The present mass model aims to characterize the 
present-day description of the distribution of matter in our Galaxy, which is needed to frame current crucial issues of Cosmology, 
Astrophysics and Elementary Particles. 



< 

O 



> 

(N 

in 

o 



% 



1 Introduction 

The evidence for the phenomenon explained in term of 
Dark Matter in Galaxies is outstanding. In particular, spiral 
rotation curves (e.g. Rubin et al. 1980; Bosma 1981) have 
revealed the presence in these objects of a dark "mass com- 
ponent". The current framework includes a central bulge, 
a stellar disk and an extended gaseous disk, all embedded 
in a spherical halo made by dark particles (of yet) unknown 
nature. In external galaxies, the dark-luminous mass dis- 
tribution has been derived by means of a) hundreths of in- 
dividual rotation curves (e.g. Spano et al. (2008)) b) 1000 
coadded RCs (Salucci et al. 2007) and c) various measure- 
ments of galaxy's gravitational potential (see e.g. Salucci et 
al. (201 1)). The result is that DM halos show a shallow cen- 
tral density profile, moreover, stellar disk radii, halo core 
and viral radii, gaseous and stellar masses, turn out to be all 
related and to lead to an Universal Rotation Curve (Salucci 
et al. 2007), likely the final state of complex physical pro- 
cesses governing the formation of galaxies. 

Let us now consider the Milky Way. It is logic to ask: 
what is its distribution of dark and luminous matter? Does 
it conform to that of external galaxies? Moreover, living 
"inside" the object does help us in better disentangling the 
diflTerent mass components? Also of interest, do we know 
well enough how the DM particles are distributed around 
the Galaxy, in order to predict the signals of their annihila- 
tions to be expected in current experimental surveys? 

It is not easy to address these questions and obtain a 
proper mass model of the Galaxy. The MW is an intrinsi- 
cally complex object: observations are difl&cult to interpret 



due to our location within it. We cannot measure the MW 
circular velocity, and we must resort to more indirect kine- 
matical measurements to infer its gravitational field. In fact, 
inside the solar circle (e.g. McClure-Grifl^iths and Dickey 
2007), we measure the rotating HI disk terminal velocities 
Vt that relate to the circular velocities V{r) once we as- 
sume values for 7?©, the distance of the Sun from the Galaxy 
center and V© the the value of the circular velocity at the 
Sun's position: V{r) = Vrir) + V© r/7?©. Both quantities are 
presently known within an uncertainty of 10% (e.g. McMil- 
lan, Binney 2009 that triggers a not negligible uncertainty 
in the derived circular velocity V(r). Outside the solar cir- 
cle, the circular velocities V(r) are inferred by applying the 
Jeans equation to the kinematics of populations of tracer 
stars. Remarkably, the result strongly depends on the (un- 
known) dynamical state of tracers and it is also aflTected by 
uncertainties in their photometric and kinematical measure- 
ments (Battaglia et al. 2005; Xue et al. 2008; Brown et al. 
2009). The mass modeling method features (see Caldwell, 
Ostriker (1981), and e.g. Sofue (2009)): a) an ad hoc halo 
density profile b) a number of assumptions on the dynami- 
cal status of the tracers and c) on their density distributions. 
All this adds to the complications due to the fact that the 
MW rotation curve is quite flat and therefore it is intrinsi- 
cally difl&cult to decompose in its dark and luminous com- 
ponents (Tonini and Salucci 2004). 

However, there are today good motivations to attempt de- 
riving a robust and reliable MW mass model, i.e.: to test 
claims of Gamma-ray emissions from annihilating DM par- 
ticles in high density regions of the Galaxy, to investigate 
the Nature of the dark particles from the galactic DM den- 



sity distribution, to understand how the Galaxy and the Lo- 
cal Group formed. A novel approach is warranted since in 
the past few years, modeling techniques and measurements 
have progressed. 

First, let us note that the mass distribution in the Galaxy 
is likely to be similar to that of any other galaxy of the same 
mass/luminosity/dimensions (Salucci et al. 2007). More in 
detail, the DM halo density profile in spirals is represented 
by the URC halo profile (Donato et al. 2009). The diff'erent 
profiles adopted in several previous works, (NFW and/or 
pseudo-isothermal halo) are not supported by present day 
observations in external galaxies. In addition, in consid- 
ering the NFW halo profile one should confront with the 
standard results emerging from simulations (Klypin, et al. 
2011), a procedure not always followed in previous works. 

Second, trigonometric parallaxes and proper motions of 
sources of maser emission associated with high-mass star 
forming regions in the disk of the Galaxy have been re- 
cently available (Brunthaler et al. 2010). These measure- 
ments allows one to determine source distances and proper 
motions in a direct and geometrical way, and to obtain their 
full 3 -dimensional locations and velocity vectors leading, 
for particular line of sights, to the actual Galaxy circular ve- 
locity. This provides us with a (limited) number of true and 
direct determinations of the velocity speed on the Galactic 
disk. 

Third, a very extended and large amount of radial mo- 
tions of stars, from which to extract the dispersion velocity, 
has now become available out to 80 kpc and more (Xue et 
al. 2008; Gnedin et al. 2010; Deason 2012). 

Finally, with a profile independent method, we derived 
a most conservative value of the DM density at the Sun 
Pdm{Rq) = 0.43 ±0.11+0.11 (Salucci et al. 2010). This 
measure will be used as a cross-check of the Galaxy mod- 
eling. 

The aim of this work is to obtain, despite various un- 
certainties and difficulties still present, a robust MW mass 
model of the Galaxy, at the best of the knowledge of 2013. 
In the next section we will present the adopted mass com- 
ponents, in section 3 we describe the observational data, in 
section 4 we report the fits to these data and their results; 
in sections 5, 6 we describe the consequences for direct and 
indirect DM detection, and in 7 we discuss our conclusions. 



2 Mass Components 

In modeling the Galaxy, the main mass components are: 
the central bulge, the (thin) stellar disk and the DM halo. 
The observational uncertainties and the systematic biases 
on these components, which will be discussed below, make 
negligible (for our aims) other minor mass components 
present: the central bar, the thick and HI disks, and the stel- 
lar halo. 



Bulge. The MW central bulge region has been subject to 
a number of probes which along the years have increased 
our knowledge of its mass distribution. Nevertheless, due 
to severe extinction in the very central region in diverse 
wavelengths, there remain crucial uncertainties in the in- 
ner 1-2 kpc. In particular, the state of the art (Picaud and 
Robin 2004; Bissantz et al. 2004; Robin et al. 2012) knowl- 
edge seems to point quite clearly to a two-components 
structure: a more massive 'boxy' bulge extending up to 
~ 1.5 kpc superimposed on a minor subleading 'bar', ex- 
tending up to ~ 4 kpc. The estimates of the total mass of 
these components is however very difficult, because obser- 
vations can only probe the mass density at a minimal an- 
gle of 2° from the galactic plane, i.e. is at the outer bor- 
der of the boxy bulge distribution, thus requiring some sort 
of arbitrary modeling of the distribution in order to infer 
the total mass. Kinematical probes such as surveys of the 
rotational velocities are also limited by the large velocity 
anisotropics (see e.g. Soto et al. (2012)). As a result, the to- 
tal bulge mass has to be considered as an unknown param- 
eter, subject to a lower bound which we take conservatively 
as- 0.8-1 X 10^^ Mq. 

Due to the above uncertainties and non circularities ev- 
ident in the central region we will restrict our analysis of 
the rotation curve to radii larger than 2.5 kpc. At these dis- 
tances, the shape of the mass distribution of the (dominant) 
bulge component is irrelevant, and the bulge can be safely 
be treated as a point-mass in the center. No results in this 
paper depend on the bulge density profile, as opposed to the 
value of its mass that plays a crucial role in mass modeling. 

Disk. As in most of Spirals the MW surface density of 
the stellar disk is exponential, which for the thin disk we 
parametrize as (Freeman 1970) 

Md 



ctd 



-tIRd 



InRl 



(1) 



with Rd = 2.6 ± 0.5 kpc (Robin et al. 2008). 

The value of the thin disk mass Md is rather unknown, 
although it can be constrained to lie in the range Mo = 
5-7 X 10^0 Mo (Picaud and Robin 2004; Juric et al. 2008; 
Robin et al. 2008; Reyle 2009). At the same time the disk 
scale length estimate and uncertainty may be subject to sig- 
nificant biases (substructures, coverage, etc, see Robin et 
al. (2008)). The disk thickness is constrained to be of the 
order of few hundreths of parsecs and in this work its pres- 
ence and its variations have a negligible impact on the stel- 
lar contribution to the MW rotation curve (less than few %, 
see e.g. Salucci et al. (2010)). It is therefore justified to take 
the infinitesimal thin disk approximation. At the same time, 
the thick and HI disks, whose masses are estimated to be ~ 
one tenth of the thin disk one Robin et al. (2008); Nakan- 
ishi, Sofue (2003) can be neglected in the present analysis 
and be considered to be part of the above range of Mo- 



Dark Halo. We adopt, after Salucci et al. (2007) and (Do- 
nate et al. 2009), the Universal Rotation Curve (Burkert) 
profile. In addition, as a comparison, we consider also a 
NFW dark halo. These density profiles are parametrized 
through the density scale pn and scale radius Rh 
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where x = r/Rn. Notice that Rh has a diff'erent phys- 
ical significance in the two models: the radius at which 
d log Pnfw = -2 and the radius of the region of approxi- 
mately constant density. 

3 Observational Data 

3.1 Terminal velocities 

The circular velocities inside the solar circle Rq can be "re- 
constructed" from HI terminal velocities Vn as (r/, Vt) = 
iXi, Vti + J/ Vq, where yi = tiIRq. Notice that yi is actually 
an angular measurement, with uncertainty which is negli- 
gible for the present work. The fit consists in comparing 
Vn with the terminal velocities as predicted by the model. 
In doing so, j/, is translated into r/, thus involving Rq in 
addition to V©. 

Data from a number of works (Malhotra 1995; Clemens 
1985; Alvarez et al. 1990; McClure-Griffiths and Dickey 
2007) are reported in appendix A and are binned in intervals 
of ~ 0.5 kpc, see figure 1. The binning procedure is deli- 
cate and we stress that some care has to be used in estimat- 
ing the uncertainties later used for the fit. In fact, the errors 
on single measures are very small, even below 1 km/s, but 
data themselves are widely spread, due to diff'erences be- 
tween the values coming from diff'erent regions (e.g. at op- 
posite longitudes). As a result, the uncertainties are mainly 
data-driven, and they have to be estimated accordingly. We 
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believe that binning data separately for each measurements 
provides a better description of the observative constraints 
with respect to the usual procedure of a global increase of 
the individual errors. 

3.2 Maser velocities 

Very precise measures of position and proper motions of 
maser sources, located at diverse galactocentric distances, 
e.g. Reid et al. (2009); Sanna et al. (2012); Honma et al. 
(2012) provide us with a number of circular velocities. 
While the galactocentric radii referred to these measure- 
ments are not crucial (in view of the almost flatness of 
the MW rotation profile) their amplitude is important in 
that it can resolve the degeneracy that V© and Rq have. 
The outermost of such measurements reach the galacto- 
centric radii of ~ 13 kpc and hint to a flat rotation curve, 
V{r ^ 13kpc) ^ V(7?o). They also set the scale of the cir- 
cular velocity, otherwise uncertain between 200 km/s and 
250 km/s, to the upper end, V© ~ 250 km/s 

Thus, the two important outcomes of the maser obser- 
vations (Bovy et al. 2009; Brunthaler et al. 2010) are: i) 
the quantity cd = Vq/Rq = 30.3 ± 0.9, consistent with the 
previous determinations through the motion of Sagittarius 
A*, and ii) a direct estimate of the sun circular velocity of 
Vq ^ 239 ± 7km/s, which is also in agreement with the 
knowledge of the Sun's galactocentric radius via oj (see 
e.g. Schonrich (2012)). See also (Honma et al. 2012) for 
more recent similar measurements. 

These two results play an important role in the mass 
modeling: we will use i) as a constraint between V© and 
Rq because its precision is much higher than the uncertain- 
ties in other quantities, and use ii) as the crucial measure of 
the Sun's LSR rotation velocity. 

3.3 Stellar halo velocity dispersions 

From the measured radial velocities, the velocity disper- 
sion crj^(r) of difl'erent halo populations of tracer (bright) 
stars have been estimated from the solar neighborhood out 
to ~ 80 kpc. They provide unique constraints on the rota- 
tion curve at large radii and thus on the DM density profile. 
In particular, as we will see below, the halo scale radius Rh 
is significantly related to these measurements and to their 
uncertainties, that thus require a detailed discussion. 

Assuming virialization, each population traces the gravi- 
tational potential, and we can use the spherical Jeans equa- 
tion to link the measured velocity dispersion and the Galaxy 
gravitational potential. The tracers density is well repre- 
sented by power law fit, p/ oc r'^'^ , so that the Jeans equa- 
tion can be written as 



Figure 1 : Terminal velocities, rebinned, see appendix. 
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Figure 2: Tracers velocity dispersions: "1" Gnedin 2010 
(red squares), and "2" Xue 2008 (blue circles). 



Note that we must allow different JSt for each population of 
tracers, because although they share the same gravitational 
potential there is no argument for them to have the same 
kinematics yS/ and ct/. 

The Jeans equation (4) relates the circular velocity model 
V^(r) , the velocity dispersions cri(r) and the anisotropics 
JSi. Thus, for each given velocity model V(r) the equation 
can be integrated to find the predicted velocity dispersion 
as a function of radius cri(r). For each population sepa- 
rately, the dispersion profile so computed from the model 
can then be directly compared with data. The integrated 
cri(s) is determined up to a free constant multiplied by a 
known function, solution of the homogeneous equation, 
^Lm ^ Qj(yi-^/^i)^^^r^ Considering for instance a linear de- 
pendence of the anisotropy on the radius, yS/(r) = /3i -\- /3'.r, 
the function is simple, cr^^^ = (T^ r^-2^e"^^'^, where o" is 
the free constant, defined so that (T/(80kpc) = o"/.^ In prac- 
tice, although the constants d"/ enter the fit as additional free 
parameters of the mass model, their eff'ect is very limited, 
because the dispersion kinematics is very uncertain at large 
radii (see figure 2). In practice, the eff'ect of d"/ is just to re- 
duce even more the constraining power of the data at large 
radii; So, for simplicity and without modifying the fit re- 
sults, we set d-1,2 = 105 km/s. We checked directly that also 
the derivative JS'. of the tracers dispersion anisotropy has a 
negligible impact on the fits, due to the large uncertainties 
of the observed dispersion velocities and to the requirement 
that, at large radii, /5/(80 kpc) should not be unphysical. The 
dependence on/5J can well be mimicked by a minor shift of 
the constant anisotropy yS^, and given the uncertain situation 
with the latter, we can safely setySJ = 0. 

The procedure of integrating the Jeans equation is of 
course equivalent to reconstructing the values of V^ from 
the measures of cr^, which however would lead to increased 



(and correlated) errors due to the derivative in (4). We pre- 
fer to fit the actual kinematical data. 

We use two populations, the (mostly) BHB stars with ve- 
locities surveyed by HVS (Gnedin et al. 2010), hereafter 
"1", and the SDSS DR6 BHB stars (Xue et al. 2008) here- 
after "2". Both series of binned data are taken as recalcu- 
lated in (Gnedin et al. 2010) and are displayed in figure 2. 
For simplicity we consider spherical tracers density distri- 
butions, for which one has yi ^ 4 and 72 ~ 3.5. These 
values provide a good spherical fit to BHB density and ve- 
locity dispersions as confirmed also by a recent detailed sur- 
vey, see (Deason 2012). A very recent study using again 
BHB stars from the SDSS DR8y (Kafle et al. 2012) con- 
firms values for the radial velocity dispersions practically 
coincident with the older analysis of Xue et al. (2008) that 
we use. This new analysis also provides an estimate of the 
dispersion anisotropics at various radii. For r > 16 kpc it 

points to an average anisotropy of 0.6 ± 0.5.^ This is 

obviously an extremely useful information for determining 
the outer Galactic rotation curve. 

We will not consider in this work the tracers velocity dis- 
persions at small radii r < 25 kpc, because: i) data show a 
clear break both in the density and in the velocity disper- 
sion Deason (2012); ii) the dispersion anisotropy seems to 
vary rapidly in that region; iii) at such small radii the ef- 
fect of the updated values of V© and Rq may have to be 
reconsidered in the derivation of the dispersions from data, 
and iv) at these radii the masers information on V© already 
provides a good constraint on the rotation curve. 

From the two sets of data in figure 2 and the Jeans equa- 
tion one can immediately infer some conclusion regarding 
the Pi. First, we note that both velocity anisotropics ap- 
pear to have a fairly limited slope \d In cr^/d In r\ < 0.1, 
and their average values are cri ~ 100-110. In the limit 
of flat rotation curve and considering the recent estimate of 
Vq - 240 km/s (see previous section), one readily obtains 
from the Jeans equation (4) jSt ^ -0.5-0.8. So, the i.e. 
fairly tangential halo velocity dispersions. This agrees well 
with the evidence of a tangential outer halo (Kafle et al. 
2012). Actually (see below) the preferred anisotropy would 
be even more tangential, i.e. P2 < -1, due to the high Vq. 
In this respect, a decreasing rotation curve in place of a flat 
one at large radii can help relieving this tension. 

As a result, the value of Rh resulting from the fit will be 
correlated withjSi, with higher /?/ corresponding to smaller 
Rh' This will turn out to be the dominant source of un- 
certainty in the predicted values of Rh, as discussed below. 
The measure of 132 in (Kafle et al. 2012) is still preliminary. 
Accordingly, we will usq/32 - -0.5+0.5 as a fiducial range, 
keeping in mind possible updates in the future. 



^We prefer to state explicitly this boundary condition at finite radius, 
rather than considering vanishing cr/ at infinity, which is not observation- 
ally motivated and may well turn out to be untrue. 



^In Deason (2012), the velocity dispersion anisotropy at large radii was 
estimated to bej32 ~ +0.5, by using an incorrect lower rotation curve as a 
prior knowledge, in particular with V© ^ 220km/s. 



Moreover, we note that the dispersion velocities of pop- 
ulation 1 (red squares) lie on average above those of 2 (blue 
circles). Because at the same time 71 ~ 4 is larger than 
72 ~ 3.5, to trace the same gravitational potential i.e. to sat- 
isfy eq. (4) with a common V^(r), the only possibility is to 
havcy^i -yS2 ~ 0-5 or more. This is formally possible, but 
since the two samples are not relative to radically different 
populations, it is unlikely that their dispersion anisotropics 
differ so much. Actually, this difference may be a system- 
atic bias, probably linked to the different ways the disper- 
sions are estimated from observative data.^ Therefore, we 
prefer to leave some tension in the fit than stretching the 
anisotropy values, and adopt /^i = P2 + 0.5. 

4 Fit 

Each choice of model parameters determines V{r), Vrir) 
and by means of the integration procedure described above, 
the velocity dispersions cr^(r). 
We fit 

• the binned terminal velocities of figure 1 against those 
calculated from the modeled rotation curve . 

• the tracers velocity dispersions of figure 2 against the 
dispersion profiles calculated from V{r). 

• the circular velocities of the Maser regions. 

• Vq as determined above from maser observations. 

Denoting with xt all these different quantities (evaluated at 
their radii) and with xt^exp their observational constraint, a 
standard ;^^ = Z^iU/ - ^Uexpf I ^x^t^exp ^^ employed, and 
with our data, N = 40. 

The complete set of parameters are {p//, Rh, Mb, M^, 
Rd, Rq} which specify the modeled galaxy rotation curve, 
plus the anisotropics {y^i, yS2} of the two populations which 
determine the dispersion velocity profile. However, in prac- 
tice the number of relevant parameters can be reduced. In 
fact, 

• As discussed in section 3.3, we use/^i -132 = 0.5, and 
study yS2 in the range yS2 = -0.5 ± 0.5; 

• Rd is constrained by observations, as discussed in sec- 
tion 2; 

• Mb and Md are free parameters that lie between a min- 
imum and a maximum value. As we will see below the 
fits prefer often the lowest values for them so the min- 
imum is on the boundary, and they play effectively a 
minor role in the minimization. 



^Considering for instance the 10-20% overall change in the veloc- 
ity dispersions resulting from different estimation methods, as reported 
by Brown et al. (2009). 



• In view of the quite precise knowledge of the sun's 
angular velocity, Vq/Rq ^ 30.3 ± 0.3km/skpc as re- 
called above, we express Rq in terms of V© by means 
of the above relation. 

Finally, through the mass modeling there is a one-to-one 
correspondence between pn and the circular rotation at the 
sun Vq = y(ro) (all the other parameters kept fixed), and it 
is technically convenient to trade pn for V©, which is used 
in its place as independent variable. 

4.1 Galaxy mass modeling results 

The fits have been performed with values of /32 in a range 
of 132 € [-1,0], and we report in figure 3 the resulting best 
fit halo parameters and confidence intervals for each case 
of yS2 = 0, -0.5, -1 (left to right in each frame). In figure 4 
we show the central (preferred) best fit Galactic rotation 
curves, and in figure 5 we show the sections of ;^^, describ- 
ing also the effect of the other astrophysical parameters and 
their correlations. 

As discussed above, given the estimate of 132 - -0.5 
in Kafle et al. (2012), we consider this as our preferred fit 
point, and the other two as extremal cases. With these three 
choices, the values of the best reduced ;^^^^^4q are respec- 
tively 0.59, 0.41, 0.35 for the Burkert profile, and 0.9, 0.46, 
0.35 for the NFW profile. The fit is further improved with 
even more negative values of ft, a fact which hints to pos- 
sible tension between data and the model. Thus a careful 
reconsideration of the uncertainties in the tracers velocity 
dispersions could lead to a softening of this problem. Nev- 
ertheless, due to the low values of x])of < 1' we can say 
that both models are able to provide reasonable fits to the 
data. 

From figure 3 we note that, as anticipated, the predicted 
values of Rh are strongly dependent on the values of the 
dispersion anisotropy yS2. Its possible variation is thus the 
main source of uncertainty, driving the variation in Rh and 
the ranges presented in table 1 for the halo parameters and 
derived quantities. The observational (or theoretical) like- 
lihood profile of 132 is not known at present and thus one 
can not give a definite evaluation of confidence intervals by 
marginalizing on this dominant uncertainty. Nevertheless, 
the adopted range in [32 allows us to estimate the expected 
ranges of parameters, in a consistent picture fitting all ob- 
servational constraints, by extending the 95% C.L. fit un- 
certainties in a flat range of 132- 

In figure 4 we display the best fit rotation curves, where 
the fitted terminal velocities, velocity dispersions, and 
maser circular velocities are compared with data. In fig- 
ure 5 we display the cross correlations between the fit pa- 
rameters Md, Mb, Rh, V©, by plotting the sections of ;^^ 
through the best fit point, as well as the constraints on astro- 
physical parameters. As a result, the best fits are sometimes 
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Figure 3: Best fits and confidence regions in the ph-Rh plane at 1, 2, 3cr C.L. (black, red, orange, yellow, two-parameter 
contours) for the Burkert (left) and NFW (right) profiles. In each frame the three regions correspond to tracers dispersion 
anisotropy y^o = 0, -0.5, -1 (left to right). The superimposed blue (solid and dotted) contours mark Mv/r[ 10^^ M©] = 0.1- 
10 and the shaded blue region is disfavored by estimates of the total MW mass. The thick green dotted line shows the 
Cvir-Myir relation as resulting from cosmological simulations. Parameters other than Rh and pn are fixed at their best fit 
preferred values (see text). 



constrained to lie just outside these regions. In table 1 we 
report the halo parameters according to the preferred fit and 
also quote parameter uncertainties. 

From the parameter variations in table 1 and from fig- 
ures 3 and 5, we see that the Milky Way mass modeling 
is still quite uncertain, and this is mainly due to our igno- 
rance of the dominant baryonic mass components (bulge 
and disk) and even more of the kinematics in the outer 
galaxy. The latter may be expected to improve in the 
near future due to more extensive stellar surveys and maser 
probes, on the other hand, the knowledge of the mass en- 
closed in the central region appears to be essentially ham- 
pered by the severe extinction at the relevant wavelengths. 
Still, some messages emerge from the fit. 

The best fit Burkert model features a total halo mass of 
M^ir ~ 10^2 Mo, a central DM density of ~ 4x 10^ Mo/kpc^ 
and a core radius of ~ 9 kpc. The emerging mass model 
results very similar to that of the URC relative at spirals 
like the Galaxy, although in the two cases the gravitational 
potential is obtained in completely diff'erent ways. The 
variations of j32 can significantly aff'ect the mass model- 
ing, and drive the total mass also to values larger than 
2 X 10^^ Mq. The Burkert best fit values in table 1 give 
LogCpz/T^/z/Mopc^) ^ 2, in perfect accordance with the 
trend obtained for other galaxies in a wide span of mag- 
nitudes (Donato et al. 2009). This is also the correlation 
between pn and Rh which can be inferred from the left fig- 
ure 3, so that this relation is stable under variations of yS2- 

We also analysed the model featuring NFW DM profile 
plus baryons, shown in the right panel of figure 3, as well 
as in the lower panels of figure 5 and 4. In figure 3 we also 



show the contours of the total MW virial mass M^tr and 
the concentration parameter c^ir on the right axis (using an 
overdensity parameter Ay/r = 200). 

As one can see the allowed range of halo scale radii turns 
out to be quite wide, Rh = 10-30 kpc. In addition for NFW, 
unlike the Burkert case, the correlation between pn and Rh 
in the right panel of figure 3, is of the form pnRJj -const., 
which reflects just the velocity constraints along the curve 
(both maser velocities and halo dispersions), in practice 
just a byproduct of the parametrization of the DM density 
through ph, devoid of deep significance. Similarly for the 
Vq-Rh correlation in the central frames of lower figure 5. 

In addition, two more remarks can be made. First, the 
NFW plus baryons best fit models require minimal bulge 
and disk masses as well as the largest disk scale radius, all 
at the boundary of the constraining intervals. This can be 
clearly appreciated in figure 5. If one fixes the more rea- 
sonable value of Rd - 2.5 kpc, the NFW fit worsens to 
XJ)OF ~ '^•'^' corresponding to a probability of less than 
10""^ (the Burkert fit instead raises to Xdof ~ ^-^ which 
is still very good). Thus, the acceptable fits of the NFW 
plus baryons model are only achieved by stretching the as- 
trophysical parameters to their allowed values. 

Second, a fairly high value for Cyir - 20 ± 10 is ob- 
tained. Current ACDM cosmological N-Body simulations 
imply the relation c^tr = 9.6Mw,/(/z-^ lO^^M©) (Prada 201 1) 
which is also depicted on figure 3 for comparison. The 
Galaxy reasonably has M^ir > lO^^M©, so we should have 
found Cy/r < 10, a value significantly lower than the best 
fit one. If we impose the above Cy/r-mass relationship the 
Milky- Way mass model becomes an unacceptable represen- 
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Table 1: Best fit halo parameters and uncertainties for the 
Burkert and NFW profiles. The best fit values are relative 
to/^i = -0.5, while the reported ranges correspond to 2 cr 
intervals (95.45% C.L.) calculated in the two-parameter 
space Vq-Rh (see figure 5) by considering also the range 
ySi = [0,-1]. As a result, quoted uncertainties should not 
be considered as statistical gaussian intervals, but only as 
estimates of expected values. 



tation of data. 

It is important to stress that the uncertainty in the trac- 
ers anisotropy discussed above can not be invoked to re- 
lease this tension, since a) it would require very negative 
values such as ~ -3, i.e. the existence of very tangential 
motions at odds with results of cosmological simulations, 
and b) the total MW mass would end up being unaccept- 
ably high, ~ 10^^ M©. Considering that the NFW best fits 
are already obtained by stretching astrophysical parameters 
to their extrema, models that solve the high concentration 
issue would also be too contrived. 

For these reasons, we can conclude that present data tend 
to prefer a cored profile with respect to a cusped one. Part 
of this preference can be tracked down to the terminal ve- 
locities, which require the lowest possible mass enclosed in 
the inner region (r < 8 kpc), given that the disk and bulge 
already provide enough gravitational potential (see 5)."^ Ac- 
ceptable fits can still be achieved, but at the price of finding 
a solution in tension with the ones arising within the ACDM 
paradigm. 

The value of the DM density at the Sun's radius p© is also 
an outcome of the fit, and it turns out to be consistent with 
the profile independent determination derived in Salucci et 
al. (2010) (see discussion below). 



In fact, we performed a similar analysis for the Einasto profile, which 
yields even worse fits in the inner region. 



Recalling that the main source of uncertainty in our 
fit comes from the poor knowledge of the outer tracers 
anisotropics jSt and the new recent knowledge on the lat- 
ter, it is useful to compare our findings with the results of 
some other works. 

Catena and Ullio (2009), built MW mass models in the 
process of estimating of DM local density. They have con- 
sidered both NFW and URC Burkert dark matter halos and 
tested them with a number of dynamical observables for 
the Galaxy (similar but not coincidental with those in the 
present paper, in particular with a diff'erent likelihood func- 
tion). Both models have been found to agree with their ob- 
servational data, with values of parameters somewhat simi- 
lar but coincidental with ours. Among those an higher value 
of Cyir. In Catena and Ullio (2012), with a similar global fit, 
they study also the DM phase space distribution. In view 
of the uncertainties in the halo tracers dynamics and also 
in the other available observative constraints, we consider 
quite premature such an analysis. 

locco et al. (2011) assuming state-of-the-art models for 
the distribution of baryons in the Galaxy, used microlens- 
ing and dynamical observations of the Galaxy to show that 
these data are in good agreement with the predictions of 
ACDM Dark Matter profiles. 

Deason et al. (2012) used 2000 distant Blue Horizontal 
Branch stars at Galactocentric distances of 15 kpc < r < 
50 kpc as kinematic tracers of the Milky Way mass distribu- 
tion. Their density was modeled as an oblate, power-law 
halo embedded within the spherical power-law DM poten- 
tial. Be means of distribution function method they ob- 
tained the power-law potential exponent and the velocity 
anisotropy of the halo tracers. The resulting outer circular 
velocity profile for the Milky Way when reproduced by a 
NFW halo, led to a high concentration value (Cy/r = 20) in 
agreement with present results. 

McMillan (2011) presented a simple method for fitting a 
MW mass model to observational constraints (not too dif- 
ferent from those we used in this paper). By means of a 
Bayesian approach they took input from observational data 
and expectations from the self consistent NFW plus disc 
model. They also forced the value of c^ir to lie in a range 
compatible with N-Body simulations, i.e. between 7 and 12. 
For the values of the disc scale lengths of 3.00 ± 0.22 kpc 
Solar galactocentric distance of 8.29 ± 0.1 6 kpc and a cir- 
cular speed at the Sun of 239 ± 5km/ s; MW stellar mass of 
6.43 + 0.63 X 10^0 Mo; a virial mass of 1.26±0.24x 10^^ ^^ 
and a local dark matter density of 0.40+0.04 GeV/cm^ their 
model fits observations. 

It is also worth to stress that the Deason, McMillan or 
locco approaches are not much diff'erent from ours, but they 
did not consider cored DM profiles. These works, with 
the present one, support the view that in the MW a NFW 
plus baryon model can reproduce data, but with somewhat 
stretched values of the astrophysical parameters. 
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Figure 4: Best fit rotation curve, terminal velocities and velocity dispersions for the URC Burkert profile (upper panel) 
and for the NFW profile (lower panel), compared with the adopted observational data on terminal velocities, maser data 
and halo tracers velocity dispersions (left to right). Dotted curves show the bulge, disk and DM halo separate contributions 
(left to right) 
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Figure 5: Two dimensional sections oix^ for the URC Burkert (upper) and NFW (lower) profiles, passing through the 
best fits, with fixed ^2 - -0.5. The shadowed regions mark the values of Mb (green left vertical band), Mj) (blue 
bottom horizontal band) and Rj) (the other grey bands), disfavored by bulge and disk studies. The contours mark t^x^ - 
(2.3, 6.18, 11.8), i.e. the two-parameter confidence levels of (68.27%, 95.45%,99.73%). 
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Figure 6: Galactic escape velocity as a function of galac- 
tocentric radius r for the best fit URC Burkert (black con- 
tinuous) and NFW (blue dashed) models, with their 95% 
uncertainties in the halo parameters, including the uncer- 
tainty in the anisotropy y52. The dotted line marks the Sun's 
radius. 

5 Direct DM search: local density and es- 
cape velocity 

The derived mass models provide us with the local 
value of DM density. For the Burkert profile one finds 
pQ ^ 0.49!^;^^ GeV/cm^ for the NFW values of ^ 
0.47^QQgGeV/cm^. These values are clearly consistent 
with the profile independent estimate of Salucci et al. 
(2010), as it must be since the gauss law and formula (11) 
there are satisfied. Nevertheless, the present global MW 
mass modeling has helped in reducing the uncertainties, 
from 0.2 to roughly O.lGeV/cm^, because the global fit 
constrains the values of M^ , Rq/Rd, V© and rotation-curve 
slope more strictly than the set of the upper-lower limits 
taken in Salucci et al. (2010). The central value is also 
slightly higher, due to the increased rotation speed at the 
Sun Vq, which is now better known. 

From the mass models one can derive the Galactic es- 
cape velocity, whose value at the Sun's location is impor- 
tant for direct searches through DM nuclear recoil, because 
it aff'ects the DM velocity distribution. We plot in figure 6 
the escape velocity as a function of radius, out to 300 kpc. 
One can observe that at the Sun's location the Burkert and 
NFW best fits lead to quite similar estimates of the es- 
cape velocity, we have V^^^^(Rq) - 576 ± 124km/s, and 
y^f^(/?0) ^ 613 ± 114km/s. These values are quite large 
and reflect the uncertainty due to the tracers anisotropy yS2- 
They are however in agreement with recent estimates of the 
local escape speed using RAVE data (Smith et al. 2007). 
Considering the overall variation of escape speed among 
the two models, one finds a large uncertainty of ~ 250 km/s, 
i.e. roughly 50%. This should be taken into account in di- 
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Figure 7: The "prompt" emission factor from DM annihila- 
tion as a function of the Galactic longitude, for the best fit 
Burkert (black continuous) and NFW (blue dashed) models, 
with their 2cr regions (95.45% C.L.). Note that for gamma 
rays in for instance the Fermi detector, below < 0.1-1° 
the point-spread function would smear the observed profile, 
making it eff'ectively cored for any DM profile. 



rect search experiments, especially in the regime of light 
"WIMP" DM mass of the order of ~ 10 GeV or less, where 
one relies on the upper tail of the velocity distribution to 
estimate the number of expected recoil events (see e.g. the 
analysis in McCabe (2010), Catena and UlHo (2012)). 

The uncertainty in the escape velocity at large radii r > 
25 kpc is also considerable, (i.e. from 300 to 600 km/s at 
50 kpc) and this has an impact when estimating the halo 
dispersion velocities, through the selection of outliers ver- 
sus acceptable halo tracers. Clearly in the future a modeling 
with this feedback would give a more self-consistent picture 
(see e.g. (Brown et al. 2009) for such an attempt). 

6 Indirect DM search: annihilation 

The flux from DM annihilation is conveniently expressed 
in terms of the "prompt" emission factor 



J ann\y) 



P^Rq Jl.o.s. 



pI(x) dx , 



(5) 



which we normalize by using p© = 0.4 GeV and ^o = 
8.3 kpc. The factor /^„„, as a function of the longitude £ 
from the galactic center, traces directly the angular pro- 
file of the dominant observed flux from annihilation into 
gamma rays. For annihilation into other (charged) parti- 
cles, which are then subject to bremsstrahlung, scattering 
with ISR, and nontrivial galactic difl'usion, see Cirelli et al. 
(2011). 

In fig 7 we plot Jann for the URC Burkert -\- baryons and 
for the NFW -\- baryons models. We see that for each mass 
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model the uncertainties in the galactic parameters lead to 
variations of the expected flux of a factor of ~ 5, in the 
innermost region. On the other hand, in direction of about 
40-60° from the galactic center the flux is predicted within 
a factor of 2 only, and independently of the profile chosen. 

In fact, the diff'erences between the two diff'erent mass 
models emerge only at ^ < 15°, and with a clear dis- 
criminating power only below £ < 10°, corresponding to 
r < 1.5kpc, i.e. inside the bulge region. One should thus 
bear in mind that this plot extrapolates the DM density pro- 
file in the very central region where observations can not 
constrain it. Indeed, the DM density in the bulge region 
or at shorter scales may well deviate from purely cored or 
NFW profile, still without modifying the present global fits. 

For instance, if one is willing to consider the scenario in 
which a DM core results from baryon feedback (mainly su- 
pernovae explosions) which erases the density cusp during 
galaxy formation, the same mechanism may well leave a 
'mini-cusp' in the central region, which would give a small 
contribution to the total mass inside the solar circle, but de- 
pending on the very inner density slope, may contribute to 
an evident annihilation signal from the inner zone, with a 
very localized source region of at most few degrees in an- 
gular size. 

7 Conclusions 

In this paper we have derived a model for the DM in the 
Milky Way, by adopting state of the art inner terminal ve- 
locities, probes of maser star forming regions, and outer 
stellar halo velocity dispersions. We performed our analy- 
sis for both a cored (Burkert) and cusped (NFW) profiles. 

Our fit confirms that the presence of a dark component in 
the Galaxy is beyond any doubt. 

We derived the preferred ranges of DM halo parameters 
pH, Rh and their correlation in figures 3, and reported the 
values and derived parameters in table 1 The results show 
that the DM halo is still subject to a considerable degree of 
uncertainty. This is due to the limited constraining power 
of the available observations, especially in the outer stellar 
halo region, and in the Galactic center, bulge region. The 
main source of uncertainty has been shown to be linked to 
the unknown velocity dispersion anisotropy of the tracers 
in the stellar halo, which is still very poorly constrained by 
observations, and which drives the main uncertainty in the 
halo scale radius Rh- Nevertheless, thanks to recent more 
precise measurements of the rotation speed at the Sun's po- 
sition V© via maser emission in high star- forming regions, 
the fitted halo parameters permit to draw some interesting 
conclusions. 

First, if one considers a cored profile like the adopted 
URC-Burkert one, very good fits are obtained, which can 
accommodate all available observations. The uncertainty 
in the tracers velocity anisotropy is still driving the Halo 



parameters uncertainty, resulting in the ranges reported in 
table 1, with Rh = 9.26+^;^ kpc and p// = 4.13+f^ with 
quite a strict correlation of PhRh ~ const, even varying 
the tracers anisotropy /52. The value of phRh - 10^ is in 
excellent agreement with Dark-Matter fits in many other 
similar external galaxies, and constitute a reference point 
on which one may base further analyses of the Dark Matter 
in our Galaxy. 

We analyzed also the case of NFW profile and showed 
that, while it also can be used to fit the observations, this 
was only possible by stretching other astrophysical param- 
eters to their allowed limits; e.g., the disk scale length had 
to be set to 3.1 kpc which is at its marginal upper limit, and 
the bulge and disk masses had to be as small as possible. 
This is due to the profile of inner terminal velocities, which 
require the lowest mass as possible in the central region, 
where instead the NFW cusp is. Indeed, a similar analy- 
sis for the Einasto profile showed a further worsening of 
this tension. In other words, the present data rather prefer 
a cored DM distribution, though the presence of standard 
NFW halo cannot be excluded and, if one decides to adopt 
it, a great care must be taken in choosing the values of the 
model parameters. 

Also, an important evidence is that, especially due to 
the high V© value, the NFW halo turns out to have a high 
Cyir ~ 20, at odds with results from numerical simulations 
in the ACDM setup. Of course this might be due to the 
eff'ect of dirty baryonic physics. Recent simulations have 
started to account for the eff'ect of baryons which have a 
twofold eff'ect, namely, to: i) adiabatically contract the halo, 
increasing the inner halo slope even beyond 1-2, and ii) ex- 
pel the inner dark matter due to supernovae feedback. It is 
not clear whether the net final eff'ect is to reduce the con- 
centration cjr even further with respect to the bary on-less 
case. In this scenario a DM profile much different from the 
original NFW one is likely to emerge, and in our opinion to 
converge with the URC one. 

In the light of the fitted halo parameters, we analyzed 
also the predicted Galactic escape velocity profile and the 
flux from DM prompt annihilation. The escape speed is 
found to be very uncertain, independently from the pro- 
file chosen: at the Sun's radius it ranges in the interval 
500-750 km/s, which has an impact on direct dark matter 
searches, especially for low (WIMP) DM masses of the or- 
der of lOGeV. The escape velocity turns out to be quite 
uncertain also at large radii in the outer galaxy, with an im- 
pact on the removal of unbounded outliers from the samples 
of stellar halo tracers. 

Regarding the flux from DM annihilation, we assessed 
the discriminating power on the cored versus cusped DM 
halo hypotheses of indirect DM searches toward the galac- 
tic center. In figure 7 we displayed the expected fluxes from 
prompt DM annihilation, which shows that a discrimination 
is possible only in the region inside 10° from the Galactic 
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center, or inside 1.5 kpc. This is right in the bulge region, 
and may be subject to a considerable astrophysical back- 
ground uncertainty. Let us also stress in this respect, that 
the observational data are able to constrain the DM profile 
only at radii larger than about 2 kpc; but as the physics lead- 
ing to the inner DM density profile is still not understood, it 
may well be the case that the real profile will deviate from 
the simple core or cusp shapes. 

The DM density at the Sun's location is also an outcome 
of the mass modeling, and the values found are Ph(Rq) - 
0.49 J^f for the Burkert profile (and 0.47 1+^f for NFW). 
They are clearly in agreement with the profile-independent 
determination Salucci et al. (2010), but its uncertainties are 
reduced by a factor of two from ~ 0.2 down to roughly 
< 0.1 GeV/cm^. This is due to the fact that in the present 
global modeling the rotation curve is more constrained by 
inner and outer galactic observations and the sun's rotation 
speed is known better. 

In the present study, we tried to assess realistic ranges for 
the DM halo parameters, but still being far from claiming to 
give precise statistical significance intervals. The reason is 
due to the lack of knowledge of realistic likelihood profiles 
for many important quantities (the stellar halo anisotropy 
parameters but also the disk and bulge masses) so that any 
arbitrary choice of ad-hoc priors would translate in unre- 
alistic or unstable results. We try to avoid such unjustified 
predictions in parameters, both as a matter of principles and 
also to avoid complicating the development of the field of 
Galactic modeling as well as the one of dark matter detec- 
tion. 



Acknowledgments 

P.S. thanks the Gran Sasso Center for Astroparticle Physics 
(CFA) and INFN, Laboratori Nazionali del Gran Sasso 
(LNGS) for hospitality during the completion of this work. 

A Terminal velocities 

We describe here the adopted data relative to terminal ve- 
locities in the inner Galactic region. We use available data 
from (Malhotra 1995; Clemens 1985; Alvarez et al. 1990; 
McClure-Griffiths and Dickey 2007) versus the absolute 
value of the longitude. Since the data show systematic vari- 
ations due to local gas motion, it is important not to trust ei- 
ther the small individual errors (often smaller than 1 km/s) 
or the high number of data points. Indeed, even an infinite 
increase in the number of measures would only amount to a 
more precise map of the gas motion, and not imply a more 
precise determination of the circular velocity. As a result, 
the adopted uncertainty in the circular velocities should fol- 
low the dispersion of the data, which directly traces the 
local variations and does not become precise with an in- 



crease of data points. We perform a binning of ~ 0.5 kpc 
and consider the weighted average and the dispersion in- 
side each bin as data for the global fit. The binned aver- 
age is performed by first binning each dataset separately (to 
avoid dominance of more populated datasets) using origi- 
nal errors enlarged by 7 km/s (to avoid the dominance of 
sets with tiny experimental error). Then, inside each bin, 
we average the results from diff'erent datasets, and consider 
the weighted average and the dispersion inside each bin as 
data for the global fit. The resulting binned terminal veloc- 
ities are shown in figure 1. In figure 8 we show also the 
original data, together with the binned velocities, where for 
convenience we have reconstructed the circular velocities 
for Vq = 240 km/s and Rq = ^ kpc. 
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Figure 8: Circular velocities reconstructed from terminal 
velocities (blue, thin) and after their binning (red, thick). 
For Vq = 242 and 7?© = 8 kpc. 
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